ANALYSIS OF THE EFFECTS OF RESIDUAL STRAINS AND DEFECTS ON 
SKIN/STIFFENER DEBONDING USING DECOHESION ELEMENTS 


Carlos G. Davila 

Aerospace Engineer, Analytical and Computational Methods Branch, Senior Member, AIAA. 
NASA Langley Research Center, Hampton, VA 23681 

Pedro P. Camanho 

Assistant Professor, University of Porto 
Porto, Portugal 


Abstract 

Delamination is one of the predominant forms of 
failure in laminated composites especially when there 
is no reinforcement in the thickness direction. To 
develop composite structures that are more damage 
tolerant, it is necessary to understand how 
delamination develops and how it can affect the 
residual performance. A number of factors such as 
residual thermal strains, matrix curing shrinkage, and 
manufacturing defects affect how damage will grow in 
a composite structure. It is important to develop 
analysis methods that are computationally efficient that 
can account for all such factors. The objective of the 
current work is to apply a newly developed decohesion 
element to investigate the debond strength of 
skin/stiffener composite specimens. The process of 
initiation of delaminations and the propagation of 
delamination fronts is investigated. The numerical 
predictions are compared with published experimental 
results. 

1. INTRODUCTION 

Interlaminar cracks (delaminations), as a result of 
impact or a manufacturing defect, can cause a 
significant reduction in the compressive load-carrying 
capacity of a composite structure. The stress gradients 
that occur near geometric discontinuities such as ply 
drop-offs, stiffener terminations and flanges, promote 
delamination initiation, trigger intraply damage 
mechanisms, and cause a significant loss of structural 
integrity. The study of delamination mechanics may be 
divided into the study of delamination initiation and 
the analysis of delamination propagation. 
Delamination initiation analyses are usually based on 
stresses and use criteria such as the quadratic 


interaction of the interlaminar stresses in conjunction 
with a characteristic distance. 1,2 The characteristic 
distance is an averaging length that is a function of 
geometry and material properties, so its determination 
always requires extensive testing. 

Most analyses of delamination growth apply a 
fracture mechanics approach and evaluate energy 
release rates, G, for self-similar delamination growth. 
The G values are usually evaluated using the virtual 
crack closure technique (VCCT) proposed by Rybicki 
and Kanninen 3 . The approach is computationally 
effective since the energy release rates can be obtained 
from only one analysis. 

In the present paper, an approach is proposed that is 
well suited to progressive failure analyses where 
delaminations are present. The approach consists of 
placing interfacial decohesion elements between 
composite layers. The proposed constitutive equations 
for the interface are phenomenological mechanical 
relations between the tractions and interfacial 
separations. With increasing interfacial separation, the 
tractions across the interface reach a maximum, 
decrease, and vanish when complete decohesion 
occurs. The work of normal and tangential separation 
can be related to the critical values of energy release 
rates 4 . 

In order to predict the initiation and growth of 
delamination, an 8-node decohesion element is 
developed and implemented in the ABAQUS finite 
element code 5 . The decohesion element is used to 
model the interface between sublaminates or between 
two bonded components. The material response built 
into the element represents damage using a cohesive 
zone ahead of the crack tip to predict delamination 
growth. 
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2. ELEMENT FORMULATION 


2.1 Kinematics 

A zero-thickness decohesion element with 8-nodes is 
proposed to simulate the resin-rich layer connecting 
two laminae of a composite laminate. The constitutive 
equation of zero-thickness decohesion elements is 
established in terms of relative displacements and 
tractions across the interface. The vector defining the 
relative displacement in global coordinates is: 

\ =u t~ u i = N k u fa- N k u ki - N k u ki [ 1 ] 

where N k are standard Lagrangian shape functions, and 
w + and uj are the components of the displacements of 
the top and bottom surfaces, respectively. 

For a general element shape and alignment, the 
normal and tangential relative displacements 8 S must 
be determined in local coordinates. The transformation 
between the global and the local coordinate fame is 
accomplished using the rotation tensor 6 si defined in 

Ref. 6: 


8 V — <9 CJ A, — O vi N k u ki — B'frUfr 


J sik u ki 


[ 2 ] 


2.2 Constitutive Equation 

The need for an appropriate constitutive equation in 
the formulation of the decohesion element is fun- 
damental for an accurate simulation of the interlaminar 
cracking process. A constitutive equation is used to 
relate the traction, z, to the relative displacement, 8 [ at 
the interface. Some softening models that have been 
proposed are shown in Figure 1. One characteristic of 
all softening models is that the cohesive zone can still 
transfer load after the onset of damage in Figure 1). 
For pure mode I, II or III loading, after the interfacial 
normal or shear tractions attain their respective 
interlaminar tensile or shear strengths, the stiffnesses 
are gradually reduced to zero. The area under the 
traction-relative displacement curves is the respective 
(mode I, II or III) fracture energy. Using the definition 
of the J integral 7 , it can be shown that for small 
cohesive zones 4 : 


J* r(d)dd = G C [3] 

where G c is the critical energy release rate for a 
particular mode, and 8 f is the corresponding relative 
displacement at failure (8 f =8 pp , 8 pro , Sn n , S Ne , or 8 reg in 
Figure 1). 
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Linear softening (lin) 
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Regressive softening (reg) 

Needleman (Ne) 



Figure 1. Softening constitutive models. 


2.3 Mixed-Mode Delamination Criterion 

In structural applications of composites, 
delamination growth is likely to occur under mixed- 
mode loading. Therefore, a general formulation for 
decohesion elements must deal with mixed-mode 
delamination growth problems. 

Under pure mode I, II or III loading, the onset of 
damage at the interface can be determined simply by 
comparing the tractions with their respective 
allowable. Under mixed-mode loading, however, 
damage onset may occur before any of the stress 
components involved reach their respective allowable. 
The mixed-mode criterion proposed assumes that 
damage initiation can be predicted using the quadratic 
failure criterion: 





= 1 


[4] 


where z 3 is the normal traction, and z 2 and Ty are the 
transverse tractions. T and S are the nominal normal 
tensile and shear strengths, respectively. The operator 
<x> is defined as x if x>0 , and 0 otherwise. 

The delamination mechanisms in mode II and mode 
III are assumed to be the same. Therefore, mode III 
can be combined with mode II by using a total tan- 
gential displacement 8 skea r defined as the norm of the 
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two orthogonal tangential relative displacements Si 
and S 2 . 


Sshear-iW+Sl [ 5 ] 

The total mixed-mode relative displacement S m is 
defined as: 


^ 3 ) +$ shear 


[ 6 ] 


The exponent a in the power law is usually selected 
to be either 1 or 2 , in which case the criterion is a two- 
parameter interaction law with parameters G IC and 
G IIC . A recently proposed criterion, the B-K criterion 8 , 
is established in terms of the single-mode fracture 
toughnesses G IC and G IIC and a parameter 77 : 


G T - G Ic + ( G IIc~ G Ic) 


V u r; 


[ 12 ] 


where S 3 is the relative opening (mode I) 
displacement. Using the same penalty stiffness K P in 
modes I and II, the tractions before the onset of 
damage are: 


r i =K P S i ,i = 1,2,3 


[7] 


The single-mode failure initiation displacements are 
then: 


cO _ cO _ oO _ *5 
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[ 8 ] 
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If the relative opening displacement S 3 is not zero, 
the mode mixity can be expressed by: 


p= 


shear 

$3 


[9] 


The mixed-mode damage initiation displacement is 
obtained by substituting Eqs. 5-9 into 4, which gives: 


cO _ cO oO 

v m S3 
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[ 10 ] 


The criteria used to predict delamination 
propagation under mixed-mode loading conditions are 
generally established in terms of the energy release 
rates and fracture toughness. The most widely used 
criteria to predict the interaction of the energy release 
rates in mixed-mode is the power law given by the 
expression: 
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Reeder 9 applied the B-K criterion to mixed-mode 
test results of AS4/3501-6 and obtained a best fit using 
77=1.45, which is shown in Figure 2. 

The linear (a=l) and quadratic (a=2) criteria, which 
are also shown in Figure 2, do not correlate well with 
the experimental results for this material. Therefore, in 
order to accurately account for the variation of fracture 
toughness as a function of mode ratio, the B-K 
criterion is implemented here. 



Figure 2. G c versus G n /G T mode ratio for AS4/3501-6. 


For the bilinear traction-relative displacement 
softening law assumed here, the critical energy release 
rates in mode I and mode II are: 
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[13] 
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where S( and are the ultimate opening and 
tangential displacements, respectively. The mode I and 


mode II energies released at failure are computed 
from: 



where S^J and 8^/ are the mode I and mode II rela- 
tive displacements at failure under mixed-mode 
loading. Using the definition of the mixed-mode rela- 
tive displacement 8 m in Eq. 6 and the mode mixity 
ratio /3 given by Eq. 9, one can solve for the ultimate 
relative displacement 8^ . For the B-K criterion in Eq. 
12, the ultimate relative displacement is then: 


Si =- 




Gic + (fine ~ Gjc ) 


P 2 } 

77 




[15] 


Having defined the relative displacement corre- 
sponding to damage onset, 8^ , and total decohesion, 
8'Jn , the constitutive equation can be defined as: 


U D sr 8 r 


[16] 


The constitutive operator D sr is defined in Ref 6 as: 


D „ = 


8 sr K p , 8r<8° m 

(1 -d)K p +dK p S si ^- 
— o-> 


8 sr K p 
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[17] 

where 8 sr is the Kronecker delta, 8™ x is the 
maximum value (over time) of the relative 
displacement, and 


d= ^M> — 

cmax I c f cO 
°m \0 J m ~ O m 


max _ cO j 

m m K rfe[0,l] 


[18] 


The element stiffness matrix is obtained using the 
Principle of Virtual Work: 


^ kizv ^zv fki 


[19] 


with: 


Kkizv ~ \ A D sr B n 


'dB 


spy 


duu 


U yp + B sik 


dA 


[ 20 ] 


The irreversibility of the damage process is taken 
into account by using the maximum value of the rela- 
tive displacement, 8™ x , rather than the current value. 

The integration is performed numerically using a 
Newton-Cotes integration, which has been shown to 
perform better than Gaussian integration in problems 
involving strain softening. 4 The element proposed is 
implemented in ABAQUS 5 as a user-defined element. 


3. ELEMENT VALIDATION 

Before simulating structural components, the 
proposed approach is validated by comparing the 
predictions with experimental data obtained for 
double-cantilever beam (DCB-mode I loading), end- 
notched flexure (ENF-mode II loading) and mixed- 
mode bending (MMB-mixed-mode I and II) tests in 
unidirectional AS4/PEEK carbon-fiber reinforced 
composite. The geometry of the specimens and the 
material properties are presented in Ref. 6. The 
experimental tests were performed at different 
G I ]/(G I +G]j) ratios, ranging from pure mode I loading 
to pure mode II loading. 

The initial size of the delamination is simulated by 
placing open decohesion elements along the length 
corresponding to the initial delamination of each 
specimen. These elements are capable of dealing with 
the contact conditions occurring for mode II or mixed- 
mode I and II loading, therefore avoiding 
interpenetration of the delamination faces. The 
numerical predictions and the experimental data for the 
test cases simulated are shown in Fig. 3. 

It can be concluded that a good agreement between 
the numerical predictions and the experimental results 
is obtained. The largest difference (8.3%) corresponds 
to the case of an MMB test specimen with G II /(G I +G I] ) 
= 20 %. 
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Figure 3. Experimental and predicted load- 
displacement relations. 
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Figure 4. Skin-stiffener specimen configuration (not 
to scale). 


4. ANALYSIS OF SKIN-STIFFENER 
DEBQNDING 

Most composite components in aerospace structures 
are made of panels with co-cured or adhesively bonded 
frames and stiffeners. Testing of stiffened panels has 
shown that bond failure at the tip of the stiffener flange 
is a common failure mode. Comparatively simple 
specimens consisting of a stringer flange bonded onto 
a skin have been developed by Krueger et al. to study 
skin/stiffener debonding 10 . The configuration of the 
specimens studied in Ref. 10 is shown in Figure 4. The 
specimens are 203 mm-long, 25.4 mm- wide. Both skin 
and flange were made from IM6/3501-6 
graphite/epoxy prepreg tape with a nominal ply 
thickness of 0.188 mm. The skin lay-up consisting of 
14 plies was [0/45/90/-45/45/-45/0] s and the flange 
lay-up consisting of 10 plies was [45/90/-45/0/90] s . 

The properties of the unidirectional graphite/epoxy 
and the properties of the interface reported in Ref. 10 
are shown in Tables 1 and 2, respectively. 

The relative mode I and mode II displacements for 
damage onset are 8 ® =61xl0" 6 mm and 
82 =68xl0" 6 mm, respectively. The corresponding 
ultimate relative displacements are 

<5/=2.46xl0" 3 mm and £/=17.5xl0" 3 mm. The 

parameter 77 = 1 .45 for the B-K criterion is taken from 
test data for AS4/3501-6. 


Table 1. Material Properties of IM6/3501-6 
Unidirectional Graphite/Epoxy 10 


E 11 

(GPa) 

E22 - E 33 V12-V13 

(GPa) 

v 23 

Gl2“G 13 

(GPa) 

G23 

(GPa) 

144.7 

9.65 0.3 

0.45 

5.2 

3.4 

Table 2. 

Properties of the interface 10 


G Ic (N/mm) G„ c (N/mm) T (MPa) 

S (MPa) 

h 

0.075 

0.547 

61 

68 

1.45 


A complete analysis of the delamination growth in 
the tension specimen requires a high degree of 
complexity. For instance, Krueger 10 developed highly 
detailed two-dimensional models using up to four 
elements per ply thickness. The approach taken here, 
on the other hand, is to determine if it is possible to 
predict the debond load of the specimen using 
decohesion elements in a much coarser three- 
dimensional model. To keep the modeling difficulties 
low and the approach applicable to larger problems, 
the model that was developed uses only two brick 
elements through the thickness of the skin, and another 
two through the flange. The complete model consists 
of 1,002 three-dimensional 8 -node C3D8I elements 
and 15,212 degrees of freedom. To prevent 
delamination from occurring at both ends of the flange 
simultaneously, model symmetry was reduced by 
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modeling the tapered end of the flange with a refined 
mesh on one side and a coarser mesh on the other. This 
model does not contain any pre-existing delaminations. 

Residual thermal effects are simulated by perform- 
ing a thermal analysis step before the mechanical loads 
are applied. The same coefficients of thermal 
expansion (a77=-2.4xlO" 8 /°C and a 2 2=3 .7 xlO" 5 /°C) are 
applied to the skin and the flange, and the temperature 
difference between the curing and room temperature is 
AT=157 °C. The flange has more 90° plies than 0° 
plies and the skin is quasi-orthotropic, so it is expected 
that residual thermal stresses are present at their 
interface at room temperature. 

4. 1 Tensile loading 

The specimen configuration corresponding to tensile 
loading is shown in Figure 5. The mechanism for 
debonding under axial load is primarily by shear lag. 

Deformed plots of the finite element model imme- 
diately before and after flange separation are shown in 
Figure 5. It can be observed that only the refined end 
of the flange separates. A detailed analysis of the 
results indicated that the coarse end of the flange also 
softens, but that the separation of the flange at the re- 
fined end relieves the stresses at the coarse end. It is 
interesting to note that static tests exhibit debond of 
only one end while fatigue tests induce debonding at 
both ends 10 . 

Note that the debond growth is not symmetric across 
the width: the debond initiates on the left corner of the 
flange, as shown in Figure 5, due to the unsymmetry 
introduced by the terminated plies at the flange tapered 
ends. 

The load-extensometer measurement predictions are 
compared with the experimental data in Figure 6. The 
initiation of delamination is marked by the sharp 
breaks in the extensometer readings. Krueger et al. 10 
defined the damage initiation load as the load corre- 
sponding to the first load-drop before flange debond- 
ing. The numerical predictions and experimental re- 
sults are shown in Table 3. 


Extensometer 



Figure 5. Deformed plots of skin/flange debond model. 

It can be observed that good accuracy in the pre- 
diction of debond loads is obtained in the 2 cases 
investigated. Without thermal loads (AT=0 °C), the 
debond load is approximately equal to the maximum 
value obtained in the tests. The thermal loads (AT=157 
°C) cause an 8.5% reduction in the predicted debond 
load, which is now well within the range of the 
experimental results. The predicted stiffness of the 
specimen is also in good agreement with the 
experimental data until a relative displacement of 
approximately 0.06 mm is reached. For displacements 
higher than 0.06 mm, the experimental data shows a 
stiffening effect. This effect is due to the fact that the 
relative displacement is measured using an 
extensometer placed between two points of the 
specimen. With increasing loads, the specimen’s 
curvature also increases due to the higher bending 
moments, and the extensometer no longer measures the 
relative displacements in the global x direction shown 
in Figure 5. 


Table 3. Experimental and numerical results (tension). 



Exper. 

Analysis AT=0 °C Analysis AT= 

=157 °C 


Load 

Load 

Error 

Load 

Error 

Damage 

initiation 

20.9 kN 

NA 

NA 

NA 

NA 

load 






Flange 

debond 

22.7 kN 

23.6 kN 

4.0 % 

21.6 kN 

-4.8 % 

load 









6 


Reaction Force, kN 



Figure 6. Predicted and experimental extensometer 
measurements. 

4.2 Three-point bending 

The three-point bending test specimen simulated is 
shown in Figure 4. The mechanism for the initiation of 
debonding under three-point bending is a combination 
of peel and interlaminar shear. Deformed plots of the 
finite element model simulating the three-point 
bending load case is shown in Figure 7. Debonding 
between the skin and the flange is considered to occur 
when a through the width crack is predicted. The pre- 
dicted and experimental load displacement relation is 
shown in Figure 8. A comparison between the pre- 
dicted and experimental damage initiation and debond 
loads is shown in Table 4. 

Good accuracy in the prediction of both damage 
initiation and debond loads is obtained only when re- 
sidual thermal stresses are taken into account. Without 
thermal stresses (AT=0° C), the debond load is 
overpredicted by 22.7%. After debonding, the analysis 
including residual thermal stresses predicts a stable 
crack growth until a load of 572 N is attained. In 
contrast, the test specimens exhibited an unstable crack 
growth of approximately half of the span of the flange. 
In order to capture the unstable crack growth a more 
refined mesh, able to predict the condition of increase 
of energy release rate with increasing debond length at 
failure (neglecting R-curve effects), ^ > q , would 

be required. The stiffness of the skin-stiffener 
specimens is predicted accurately for the two 
conditions simulated. 


Table 4 Experimental and numerical results 
(3 -point bending). 


Exper. Analysis AT=0 °C Analysis AT=157 °C 



Load 

Load 

Error 

Load 

Error 

Damage 

initiation 

428 N 

517 N 

19.4% 

419 N 

-2.1 % 

load 






Flange 

debond 

463 N 

568 N 

22.7 % 

514 N 

11.0% 


load 


Deformed d=6.5 mm. 



Figure 7. Deformed plots of skin/flange debond 
model (3 -point bending). 



Figure 8. Predicted and experimental load- 
displacement relation. 


4.3 Effects of Bond Defects 

The manufacturing processes of laminated com- 
posite materials can lead to the existence of defects in 
the bonds. The debonds can be caused by air entrapped 
between plies or by inclusions in the prepreg. 
Furthermore, accidents such as tool drops, or other 
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low-velocity impact loads, often lead to debonds that 
are difficult to detect. Under the previous circum- 
stances, cracks will be present in the material before 
the application of the design loads. For damage toler- 
ance considerations, it is therefore important to assess 
the structural behavior of composite materials in the 
presence of pre-existing debonds. 



Figure 9. Location of defects. 


The procedure previously proposed is used to 
simulate the response of tension-loaded skin-stiffener 
composite structures with bond defects. The three 
cases shown in Figure 9 are investigated. The three 
cases correspond to 8.5 mm x 0.67 mm debonds 
located at different positions along the skin-flange 
interface. The defects are simulated by placing open 
decohesion elements at the corresponding positions. 
The analyses include the effects of residual thermal 
stresses. 

The predicted load-displacement relations for an 
undamaged specimen and for the three cases under 
investigation is shown in Figure 10. The predicted 
debond loads and the difference between the debond 


loads for specimens with and without pre-existing de- 
fects are shown in Table 5. 


Table 5. Predicted debond loads (tension). 


No pre-crack Case A Case B Case C 


Flange debond 
load (kN) 
Difference (%) 


21.6 


20.7 

-4.2 


17.3 

-19.9 


19.1 

- 11.6 



Figure 10. Predicted load-extensometer measurement 
for undamaged and damaged specimens 
under tension. 


From the observation of Figure 10 and Table 5, it is 
clear that both the stiffness and the debond loads of the 
specimens are reduced by the pre-existence of bond 
defects. The effects of the pre-existing defects are 
more pronounced when the defects are located at the 
specimen’s edges (cases B and C). The largest reduc- 
tion of stiffness and debond load (19.9%) corresponds 
to case B, when the defect is placed at the position 
corresponding to the onset of debonding in the un- 
damaged case. 

5. CONCLUDING REMARKS 

The debond strengths of skin/stiffener specimens 
loaded in tension and in three-point bending are in- 
vestigated using a cohesive-zone model. The good 
agreement between the experimental and predicted 
location for crack initiation indicate that the proposed 
approach can capture the onset of damage in structures 
where the exact location of an initial crack may be 
difficult to determine a priori. 


For a skin-stiffener specimen under tension, good 
agreement between the experimental results and the 
numerical predictions is obtained. The stiffness of the 
specimen, the location of crack initiation, and debond 
loads are in good agreement with published experi- 
mental data. Neglecting residual thermal stresses leads 
to a 4.8% overprediction of the debond load. 

For a skin/stiffener specimen under three-point 
bending, good agreement between the experimental 
and the numerical damage initiation and debond loads 
is obtained only when residual thermal stresses are 
taken into account in the analysis. This result indicates 
that the residual thermal stresses of the flange and the 
skin can cause a substantial reduction in the debond 
strength of the specimens when loaded in three-point 
bending. The simulation of the unstable crack growth 
following debonding under three-point bending re- 
quires a more refined mesh. 

The effects of pre-existing bond defects in a 
specimen loaded in tension are investigated for a 
8.46 mm x 0.7 mm defect located at three different po- 
sitions. In all cases a reduction of both stiffness and 
debond load are predicted. The importance of the lo- 
cation of the defect is reflected in the range of pre- 
dicted reduction of debond load, from 4.2% up to 
19.9%. It is concluded that the approach proposed here 
can be used in the failure prediction of composite skin- 
stiffener structural components when debonding is the 
leading failure mechanism. 
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